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INTRODUCTION 


Variational  principles  apply  mostly  to  boundary  problems  where  eigen¬ 
values  are  sought.  It  is  seldom  used  for  initial  value  problems  alone  where 
the  far  end  conditions  are  neither  known  nor  specified.  If  we  use  discrete 
methods  to  solve  an  initial  value  problem,  such  as  finite  difference  method, 
only  the  initial  conditions  should  be  given.  In  the  same  way,  if  we  employ 
variational  method  with  spline  functions,  we  should  not  be  concerned  with  the 
far  end  conditions.  This  paper  gives  a  procedure  to  find  a  recursive  solution 
of  an  Initial  value  problem  by  variational  methods  using  the  cubic  hermite 
polynomial  spline  functions. 

Let  us  consider  a  dynamical  system  governed  by  the  following  equation: 

L(t)ya(t)  -  -Q(t)  (1) 

with  appropriate  boundary  conditions.  In  the  above  equation  L  is  a  linear 
operator,  ya  is  the  dependent  variable,  Q  is  a  forcing  function,  and  t  is  the 
Independent  variable. 

Some  Integral  property  in  the  form  of  a  linear  functional  of  the 
variable,*  such  as  the  inner  product  of  an  adjoint  forcing  function  Q  and  the 
solution  of  Eq.  (1)  can  be  used  for  estimation. 


Glyal 


/  Qya<*t 
co 


(2) 


*W.  M.  Stacy,  Jr.,  "Variational  Methods  in  Nuclear  Reactor  Physics,"  Academic 
Press,  1974,  p.  7. 
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The  estimate  y  which  differs  from  the  solution  ya  of  Eq.  (1)  by  an  increment 
6y  can  be  written  as 

6y  “  y  -  7a  C 3) 

Then  the  estimate  y  becomes 


tb  -  tfc  -  tfc  - 

Gty]  ■  J  Qydt  -  J  Qyadt  +  /  Qfiydt 
tQ  t0  tQ 


-  G(yaJ  +  / 


tb  - 
Q6ydt 
to 


which  is  in  error  to  first  order  in  6y  and  Q. 


(4) 


THE  VARIATION  PRINCIPLE 

A  more  accurate  estimate  can  be  made  by  constructing  a  variational 
principle*  for  Eq.  (2).  By  using  the  adjoint  variable  y  as  a  Lagrange 
multiplier  for  Eq.  (1)  added  to  G[y]  we  have 

tb  - 

Jly.y]  ■  G[y]  +  /  y  (Q+-Ly)dt 
to 


tb  -  tb  -  tb  - 

-  /  Qydt  +  J  yQdt  +  /  yLydt  (5) 

to  tQ  tQ 

In  order  that  J  be  a  variational  principle  for  G  the  following  requirements 
must  be  satisfied. 

(a)  J  is  stationary  about  the  function  ys  which  satisfies  the  relation 
in  Eq .  ( 1 ) . 


L(t)ya  -  -Q(t) 


(6) 


*W.  M.  Stacy,  Jr.,  "Variational  Methods  in  Nuclear  Reactor  Physics,"  Academic 
Press,  1974,  p.  7. 


(b)  The  stationary  value  of  J  deduced  from  Eqs.  (2)  through  (5)  is 

Jly.yl  -  G[y3l  +  G[ya]  (7) 


Consider  first  the  stationarity  of  J  by  taking  the  variation 

tb  -  tb  -  tb  - 

SJ  ■  «{/  Qydt  +  /  yQdt  +  /  yLydt} 


tb  -  tb  - 

tj  m  J  5y(Ly+Q)dt  +  /  [Qdy  +  yL6y]dt 


(3) 


We  will  make  an  effort  later  to  impose  certain  conditions  in  order  that 
the  following  equality  holds: 


tb  _  tb  — 

J  yLSydt  *  /  dyLydt 


(9) 


where  L(t)  is  an  adjoint  operator. 

By  combining  Eqs.  (8)  and  (9)  one  obtains 

tb  -  tb  —  - 

5J  «  /  dy(Ly+Q)dt  +■  /  Sy[Ly  +  Qjdt  •  0 


(10) 


Since  the  variations  6 y  and  5y  are  arbitrary  it  leads  to  the  requirement  that 
the  stationary  values  y3  and  ys  must  satisfy 

Lys  “  "Q  (11) 

Lys  -  -Q  .  (12) 

Since  Eq.  (11)  is  the  same  as  Eq.  (6),  therefore  J  is  stationary  about  the 
function  yg.  Equation  (12)  is  the  adjoint  equation  in  terms  of  the  adjoint 
operator  L,  the  adjoint  variable  y,  and  the  adjoint  forcing  function  Q. 


Using  the  relation  in  Eq.  (11)  for  the  stationary  value  of  J  from  Eq.  (5) 


we  have 


tb  -  tb  _ 

Jlys  Ysl  *  /  0Ysdt  +  /  Ys^Q+LysJdt  -  G[ys]  (13) 

co  co 


Since  J  is  stationary  and  5J  ♦  0  then 


G(yg]  >  G[yaJ  (14) 

which  is  the  requirement  given  in  Eq.  (7). 

It  is  noted  that  Eq.  (10)  contains  no  boundary  terms  to  be  satisfied. 
This  bears  an  Important  point  in  the  future  discussion  of  the  initial  value 
problems . 


BILINEAR  CONCOMITANT 

The  assumed  equality  in  Eq.  (9)  is  discussed  here  by  considering  the 
following  bilinear  concomitant:1 

tb  -  tb  — 

D  -  J  yLydt  -  J  yLydt  (15) 

fco  tQ 

The  above  expression  can  also  be  written  in  terms  of  boundary  conditions 
at  t  •  t0  and  t  »  tb.  It  is  assumed  that  these  boundary  conditions  are 
assigned  in  such  a  way  that  the  above  bilinear  concomitant  is  identically 
zero,  i.e., 

D  =  0  (16) 


1W.  M.  Stacy,  Jr.,  "Variational  Methods  in  Nuclear  Reactor  Physics,"  Academic 
Press,  1974,  p.  7. 


4 


Then  Che  firsc  variations  of  D  also  vanish. 

5D  *  <5D(<5y)  +  6D(  5y)  •  0 

Since  6y  and  6y  are  independent  of  each  other,  then 

tb  -  cb  - 

SD(Sy)  »  /  6yLydt  -  /  yLSydt 


and 


tb  -  tb  — 

SD(6y)  -  /  yLSydt  -  /  <5yLydt 


(17.) 


(18) 


(19) 


Equation  (19)  is  identical  to  Eq.  (9),  which  is  the  assumed  equality 
previously.  This  implies  that  if  Eq.  (16)  is  true  then  Eq.  (9)  or  (19)  is 
automatically  true. 


INTEGRAL  OF  BILINEAR  EXPRESSION 

The  Integral  of  a  function  is  given  as 

tb 

I  *  /  i|/(yy)dt  (20) 

co 

where  \p( yy)  is  an  arbitrary  bilinear  expression^  in  the  form 

<Kyy)  *  oy'y'  +  8y'y  +  yyy'  +  eyy  (21) 


The  prime  (')  in  the  above  expression  denotes  (d/dt). 

Equation  (20)  can  be  integrated  by  parts.  Two  different  forms  of 
integration  and  end  conditions  may  be  obtained  as  follows. 


I 


"  I 


tb  -  -  tb 

yLydt  +  (ayr+Yy)y| 
to  tQ 


(22) 


2r.  Courant  and  D.  Hilbert,  "Methods  of  Mathematical  Physics,  Vol.  I," 
Interscience  Publishers  Inc.,  1953,  p.  278. 
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or 

tjj  —  cb 

I  -  -  /  yLydt  +  (ay'  +  By)y|  (23) 

to  to 

where  the  differential  expressions  are 

Ly  -  (ay' )'  -  By'  +  (ay)’  -  ey  (24) 

Ly  -  (ay*)'  +  (By)'  -  yy*  -  ey  (25) 

The  bilinear  concomitant  given  in  Eq.  (15)  can  now  be  expressed  in  terms 

of  the  function  values  and  their  derivatives  at  the  end  points  by  equating 

Eqs.  (22)  and  (23). 

-  -  -  cb 

D«  [a(y'y-y’y)  -  (y-e)yy]  i  (26) 

co 

END  CONDITIONS  FOR  THE  ADJOINT  SYSTEM 

In  order  to  satisfy  the  expression  D  =  0  in  Eqs.  (15)  and  (16)  the  end 
terms  in  Eq.  (26)  must  vanish.  Thus  it  requires 

ab(yb'yb-yb'yb)  -  «o(yo' yo-yo' yo)  -  (Yb-Sb)ybyb  +  (Yo-9o>y0yo  =  o  (27) 
Equation  (27)  can  be  satisfied  identically  If  the  end  conditions  of  the 
adjoint  system  are  proportional  to  the  end  conditions  of  the  original  system 
as  follows: 


yb 

*  (To- Bo)ky0 

(28a) 

yo 

*  (Yb"&b)kyb 

(28b) 

yb' 

*  -ab_1«o(yb-3b)kyo' 

(28c) 

yb' 

*  -a0“lab( Yo“&0)kyb' 

(28d) 

where  k  is  a  constant. 
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The  above  expressions  give  the  required  end  conditions  for  the  adjoint 
system  in  terms  of  that  of  the  original  system.  Thus  from  Eqs.  (15)  and  (16): 

cb  -  cb  — 

D  *  /  yLydt  -  /  yLydt  =  0  (29) 

co  t0 

To  summarize,  if  one  can  make  the  end  conditions  of  the  adjoint  system 
satisfy  the  relationship  in  Eq.  (28),  the  bilinear  concomitent  D  vanishes. 

The  variation  in  Eq.  (10)  is  then  valid. 

It  is  also  noted  that  the  variation  in  Eq.  (10)  has  no  far  end  terms 
which  simplify  the  computation.  This  is  because  the  far  end  terras  may  cause 
certain  difficulties  in  many  computational  schemes  on  a  number  of  variational 
methods . 

THE  FIRST  VARIATION 

Since  the  variations  6y  and  5y  are  independent  to  each  other,  we  take  the 
first  half  of  Eq.  (10)  as 

cb  -  cb  - 

5J(6y)  »  /  SyLydt  +  /  SyQdt  =*  0  (30) 

to  t0 

Equation  (30)  is  not  in  a  ready  form  for  estimation.  We  prefer  to  use  61 
which  can  be  obtained  from  the  bilinear  expression  I  given  in  Eqs.  (20)  and 
(21).  Let 

51  -  5l(5y)  +  5l(5y)  (31) 

The  first  part  of  the  above  expression  can  be  derived  from  Eqs.  (20)  and  (21) 
as 

tb 

5I(5y)  «  /  l(oy'+Yy)5y'  +  ( 0y'+ey) 5y]dt  (32) 

t0 
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Integrating  by  parts  one  obtains 


_  tb  tb  - 

Sl(Sy)  -  (ay'+yy)6y |  -  /  6y{(ay’+Yy)’  -  (By'+ey)]dt  (33) 


It  is  recognized  that  the  integrand  in  the  last  term  of  the  above  formulae  Is 
Ly.  Solving  for  the  last  term  ve  have 


cb  -  - , cb 

/  6yLydt  »  (ay'+Yy)$y|  -  <$I(oy) 


(34) 


Substituting  Eq.  (32)  into  (34)  and  then  Eq.  (34)  into  (30)  one  obtains 
6 J(6y)  -  (abyb'+YbyoJ^Vb  “  <  Ooyo’+Wo)  5yo 
,cb 


-  /  [(ay’+yy) 5y'  +  ( By’+ey) 5y)dt 

co 

cb  - 

+  J  6yQdt  -  0 
to 


(35) 


The  above  equation  contains  only  6y  and  6y*  and  none  of  the  variation  of  the 
higher  derivative  such  as  fiy"  for  a  second  order  system.  The  dependent 
variable  also  contains  only  y  and  y'  and  none  of  the  higher  derivative  such  as 
y"  for  a  second  order  system. 

ADJOINT  VARIABLE  FAR  END  VALUE  FOR  INITIAL  VALUE  PROBLEMS 


For  a  second  order  system  the  initial  values  of  the  function  and  its 
first  derivative  are  given,  i.e.,  y0  and  y0'  are  known  in  Eq.  (28).  The  far 
end  values  for  the  adjoint  system  y^  and  y^'  are  found  from  Eqs.  (28a)  and 
(28c).  Since  the  variation  of  a  constant  is  zero,  then 

syQ  ”  «yo'  ”  o  (36) 


and 


sYb  “  Syb'  ”  0 


(37) 
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The  conclusion  6yb  •  0  in  Eq.  (37)  is  important  in  that  the  first  term  at  the 
right  side  of  Eq.  (35)  vanishes.  Thus  the  coefficient  of  is  not  necessar¬ 
ily  zero.  This  implies  that  the  function  y^  and  its  derivative  y^'  at  the  far 
end  are  not  related  as  such.  By  not  using  any  local  boundary  conditions  at 
the  far  end,  the  computation  can  start  at  the  near  end  and  carry  on  in  one 
direction. 

Thus  Eq.  (35)  is  simplified  to 

-  -  ,cb  - 

«j(6y)  -  -(Y0yo+aoyo' >«y0  +  J  6yQdt 

co 

tb 

-/  {(ey+By' )5y  +  (yy+ay' )6y' ]dt  (38) 

to 

It  is  noted  that  the  above  equation  does  not  have  boundary  terms  to  be 
satisfied  at  the  far  end  at  time  tj,.  This  is  consistent  with  the  notion  of 
"initial  value  problem"  physically. 

TRANSFORMATION  OF  COORDINATES 

The  integral  sign  in  Eq.  (38)  can  be  converted  into  a  summation  sign  if 
discrete  intervals  for  integration  are  used.  Since  the  analysis  is  an  initial 
value  problem,  without  losing  any  generality  we  may  let 


co 

*  0  and  tb  * 

1  , 

(39) 

that  is  the  independent  variable 

is  within  the 

interval 

0  <  t  <  1 

(40) 

Equation  (38)  can  be  discretized 

by  letting 

5  -  Kt  -  m+1 

(41) 

0  <  5  «  1, 

0  <  t  <  1,  m 

-  1,2,...K 

(42) 

where  K  is  the  number  of  Intervals. 
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Thus 


d£  -  Kdt  dt  -  K_1d5 

The  differential  relationship  is 

dy  dy 

—  -  K  — 
dt  d£ 


or 


where 


y’  -  Ky 


(*)-  —  () 
d£ 


(43) 


(44) 


(45) 

(46) 


Then  Eq.  (38)  becomes 
5j(5y)  -  0 


*  -  K  i  - 

*  -(Yoyo+aoKy0)6yo  +  l  I  dy(®>Qirld£ 

wl  0 

K  i  .  . 

-  I  j  [(ey(m)+6Ky(m>)6y(m>  +  {(Yy(»)+aKy(m)  }K6y(m>  ]K_1d^ 

m-1  0 


(47) 


PIECEWISE  SPLINE  FUNCTIONS 

We  may  express  the  variables  y^m^  and  y (“)(£)  in  terms  of  piecewise 
spline  function  a^($)  and  the  node  point  functions  and  y(®)  as  follows. 


y(m)(?)  .  aT(5)  y(m) 

6y(m)  „  [6Y(®>]Ta(C) 

(48) 

y (“>(?)  -  aT0;)  Y<®> 

5y(®>  -  [5Y<m>]Ta(0 

(49) 

y(n)(^)  -  aT(C)Y<®> 

• 

6y(®)  -  [SY^jTaU) 

• 

(50) 

y(a)(5)  .  aT(5)Y<®> 

6y(m)  .  [6Y<®>]TaU) 

(51) 
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(52) 


m  “  1,2, ...K 
y0  -  aT(l)Y(°) 
y0  -  It(1)y(o)  (53) 

6y0  -  6Y<°>a(l)  (54) 

If  Eqs.  (48)  through  (54)  are  substituted  into  Eq.  (47)  one  obtains 

0  -  -(5Y<0>]Ta(l)[Y0aT(l)  +  OoKaT( 1) ] Y<°> 

K  i 

+  l  [  6y(h0  ]T  K“l  /  a(£)Qd£ 
m-l  0 

K  1 

-  I  («Y(“)]T  /  a(C){eir1aT(0  +  0aT(£)]d£  y(“> 
tn-1  0 

K  -  1  . 

-  I  [5Y<®>]T  /  a(C)[YaT(0  +  aKaT(£)]di-  Y<*)  (55) 

m-1  0 

This  simplifies  to 

0  -  -[6Y<°)]Ta(l)[Y0aT(l)  +  1}  ]Y(o) 

£ 

+  l  (5Y<m)]Tq(m)  -  l  [ 6Y^m> ] T  Y^®)  (56) 

m-1  m»l 

where 

q(m)  .  k-1  J1  aU)Q(5)d£ 

0 

-  Ui<m>,  q2(m).  q3(m).  q4(n)lT  07) 

and 

p(m)  «  {a(^)[e(®)K-1aT(0  +  0U)aT<£)]  +  a(  5)lY(m)aT(5)  +  o^m)KaT( £) 1  }d£ 

-  e(<“)K~1B  +  g(m>G  +  y(“>D  +  a(“>KE  (58) 

or 

[Pij(m)]  -  eGOirMbij]  +  +  i^ldy]  +  a<“)K[ei;J]  (59) 


n 


where 


Ibijl  -  /  aU)aTU)d5 

(60) 

1  • 

( ciJ  I  •  /  a(5)aT(5)d£ 

(61) 

(dij]  *  /  a(5)aT(5)d5 

(62) 

1  #  • 

[eij]  -  J  a(5)aT(C)d? 

J  0 

(63) 

CUBIC  HERMITE  POLYNOMIAL  SPLINE 


The  cubic  Hermite  polynomial  spline  is  continuous  in  the  functional 
values  and  its  first  derivatives  across  the  nodes.  Since  we  have  no  second 
derivatives  for  a(£)  in  Eqs.  (58)  to  (63),  no  higher  order  spline  is  necessary 
for  this  problem. 

The  cubic  Hermite  polynomial  gives 


aj.(5) 

- 

1-3  ^2 

a2(5) 

m 

5-2 5 2+5 3 

a(5)  - 

a3(5) 

- 

352-2S3 

a4(5) 

m 

-5*+53 

- 

- 

whose  derivatives  are 

ai(0 

m 

-6C+652 

• 

a2(5) 

m 

1-4C+352 

a(5)  - 

• 

a3(5> 

m 

65-6C2 

a4(0 

» 

-2 5+3 5 2 

_ 

(64) 


(65) 
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It  is  obvious  from  Che  above  equations  that  the  node  point  values  are 

a(0)  -  (ai(0)  a2(0)  a3(0)  a4(0)]* 

-  I  1  0  0  0  ]T  (66a) 

a(0)  *  lai(O)  a2(0)  a3(0)  a4(0)]T 

-  [  0  1  0  0  1T  (66b) 

a(l)  -  [ai<l)  a2(l)  a3(l)  a4(l)]T 

»  [  0  0  1  0  ]T  (66c) 

a( 1)  -  [ai(l)  a2(l)  a3(l)  a4(l)]T 

-to  0  0  1  ]T  (66d) 


We  wish  to  form  a  vector  whose  components  are  taken  from  the  function  and  its 
derivative  at  the  left  node  and  then  the  same  at  the  right  node.  From  Eqs. 
(48),  (49),  and  (66)  we  have 


y(m)(0) 

- 

aT(0)Y<®> 

y(m)(0) 

- 

aT(Q)Y(®) 

y(m)(l) 

m 

aT(i)Y(m) 

y(®)(l) 

m 

aT<i)Y(®) 

10  0  0 
0  10  0 

y(®)  •«  y^®) 

0  0  10 

0  0  0  1  (67) 


If  we  define 


Y(m)  .  {Yi^®)  y2^®)  Y3(®)  Y4^m^I^ 


(68) 


Then 


Yl(ffl)  -  y(®)(0)  (69a) 
Y2(m)  -  y(m)(0)  (69b) 
Y3(m)  -  y(m)(D  (69c) 
Y4(m)  »  y(m)(l)  (69d) 


The  above  Implies  Chat  Che  same  node  point  has  been  represented  by  two 
notations  as  follows 


y(ari-l)(o)  «  y(m)(i) 

(70a) 

y(m+l)(o)  -  y(®)(l) 

(70b) 

By  expanding  Eq.  (68)  for  different  m  one  obtains 

Y(0)  -  [o  0  Y3(°>  Y4^°^]t  -  10  0  y(°)(l)  y(°)(l)]T 

(71a) 

y(D  * 

[Y^1^  y2^1}  y3(1>  Y4^1^]T  -  [y^^(o)  y^\o)  y(1)(i)  y(1)(i) 

]T  (71b) 

y(m)  .  [Yj(m)  Y2^m^  Y3^m^  Y4^m^]^  * 

[y(m)(o)  y^m\o)  y^m\i)  y(m)(i)iT 

(71c) 

Y(m+1)  _  Y2^®+^  Y3^m+^  Y4^m+^]T  * 

[y(mH)(0)  y(nrK)^0)  y(mH)^1)  y(ntH)(i))T 

(7  Id) 

Thus  we 

have 

and 

Y^(n»+1)  m  y3(®) 

(72a) 

Y2(®+1)  ■  Y4^m^ 

(72b) 

for  m  » 

0,1,2,...K. 

Similar 

to  the  above  equation  one  can  prove  from  Eqs.  (50)  and  (51) 

that  the 

adjoint  variations  are 

and 

$Yj(m+l)  ,  5Y3(m) 

(73a) 

$Y2(®H)  .  6Y4^n) 

(73b) 
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METHOD  OF  SOLUTIONS 


First  we  take  the  last  term  of  Eq.  (56)  which  is 
K 

R 3  -  -  I  6Y2(m)  «Y3<m>  5Y4<m>  ]  [  Pl/m)  ]  [Y^®)  Y2(m)  Y3<n>  Y4<m>]T 

m-1 

(74) 


Using  the  relationship  from  Eqs.  (72)  and  (73)  gives 


K 

R3  -  -  l  {[pn^Y^®4’1)  +  p12(m)Y4(m-1>  +  Pi3^m^Y3^m)  +  p14(m)Y4^m>]6Y3(m_1) 

m-1 

+  lp2l^m^Y3^n,~^^  +  P22(m^4(m“1)  +  P23(m)Y3(in)  +  P24(m)Y4(,n>  ]  SY^®"1) 

+  (pj^nOY^®--1)  +  P32^m)Y4(m-1^  +  p33(m)Y3(m)  +  p34(®)Y4(m)  ]  6Y3(m) 

+  +  P42(m)Y4<m"1)  +  p43^m)Y3<Bl)  +  p44(m>Y4(m>  ]  SY4(m> }  (75) 

R3  *  -  Ipn^Y;^0)  +  pi2^)y4^®^  +  P13^^Y3^^  +  p34^^Y4^^]  6Y3^®^ 

-  [p2l^^Y3^®)  +  p22^^Y4^^  +  P23^^Y3^^  +  p24( ^Y4( ^ ]  6Y4(0) 

K-l 

**  I  {[Pll(nrH)Y3(m)  +  pl2(nri‘1)Y4(»)  +  p13(nri-l)Y3(nri-l)  +  piA(m+l)Y4(nri-l)] 
m-1 

+  [p31(m>Y3(,B-l)  +  P32(m)Y4(m-1)  +  p33(m)Y3(m)  +  p34(m)Y4(m)]}6Y3(ra> 

K-l 

-  I  {(p21(nri-l)Y3<,n>  +  p22(,“+'1)Y4(m)  +  P23(nH-l)Y3(nrt-l)  +  p24(nH-l)Y4(nH-l)  j 
m-1 

+  [p4i(m)Y3(,B~1)  +  p42ta>Y4^®"1>  +  p43(m)Y3^m)  +  p44(m)Y4(m>]  }6Y4(m) 

-  {p3l^^Y3^~^  +  p32CK>Y4(K-D  +  p33(K)y3(K)  +  P34(K)y4(K)I6y3(K) 

-  [P41<K>Y3<K-D  +  p42Y4<K-l)  +  p43(K)y3(K)  +  p44(K)Y4(K)]fiY4<K)  (76) 


It  is  noted  here  that  the  variations  at  the  far  end  are 

SY3<k>  -  5/1,  -  0  (77) 

5Y4<k>  -  6/*,'  -  0 


(78) 


by  virtue  of  Eqs.  (36)  and  (37).  Thus  the  last  two  terms  of  Eq.  (76)  drop 
out. 

It  is  again  important  to  emphasize  here  Chat  the  computation  does  not 
contain  the  condition  placed  at  the  far  end  boundary.  The  calculation  starts 
with  the  initial  conditions  and  carries  through  in  one  direction. 

The  second  term  on  the  right  side  of  Eq.  (36)  gives 

K  _____ 

R2  -  l  [qi<®>  q2(m)  q3(m)  q4(m)  ]  [  SY^®-1)  SY^®"1)  6Y3(®>  6Y4<®>]t 

m«l 

-  q1(1^6Y3<°)  +  q^^fiY^0* 

K—l  K~1  - 

+  I  [qi^®4-1^  +  q3(®)  1 6Y3^®^  +  l  [<J2^nrt~1^  +  q4^m^]6Y4(®) 
m-1  m-1 

+  q3(K)^Y3(K)  +  q^(K)fiY4^^  (79) 

The  last  two  terms  drop  out  again  by  virtue  of  Eqs.  (77)  and  (78)  . 

The  quantity  q(®)  is  again  expressed  as 

qA(m)  „  k-1  ai(OQ(®)(OdC  *  -  1. 2. 3,4  (80) 

The  first  term  on  the  right  of  Eq.  (56)  is 

Ri  -  -[0  0  SY3(0)  6Y4(°>][0  0  1  0]T{To[0  0  10]  +  OoK[0  0  01JH0  0  Y3(0)Y4<°>]T 

-  -6Y3<°>{yoY3(0)  +  OoKY4<°)}  (81) 

Combining  all  the  above  results  and  substituting  into  Eq.  (56)  we  have 

0  -  Ri  +  R2  +  R3  (82) 
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by  virtue  of  Eqs.  (36)  and  (37).  Thus  the  last  two  terms  of  Eq.  (76)  drop 
out. 

It  is  again  important  to  emphasize  here  that  the  computation  does  not 

contain  the  condition  placed  at  the  far  end  boundary.  The  calculation  starts 

with  the  initial  conditions  and  carries  through  in  one  direction. 

The  second  term  on  the  right  side  of  Eq.  (56)  gives 

K  ______ 

R2  -  l  Ul(m)  q2(m)  q3(m)  q4(m)  1 1 SY3(®-1>  6Y4(™-0  6Y3^  6Y4^m>  )T 

m-1 

-  q^D^^0)  +  q2(1)6Y4<0> 

K-l  K-l 

+  l  +  q3(m)]5Y3^m)  +  £  Iq2^nH'1^  +  q4(m>  ]  SY^®) 

m-1  m-1 

+  q3^K)5Y3^)  +  q4(K)fiY4(^0  (79) 

The  last  two  terms  drop  out  again  by  virtue  of  Eqs.  (77)  and  (78)  . 

The  quantity  q(m)  is  again  expressed  as 

qA(m)  .  K-l  f*  ajlU)Q(m)(5)d5  t  -  1,2, 3,4  (80) 

The  first  term  on  the  right  of  Eq.  (56)  is 

-  -[0  0  6Y3(0)  6Y4<0)][0  0  1  0]t(to10  010)+  OqKIO  0  01)}[0  0  Y3(0)Y4(°>]T 

-  -«Y3<0>{yoY3<°>  +  aoKY4^°)}  (81) 

Combining  all  the  above  results  and  substituting  into  Eq.  (56)  we  have 

0  -  Ri  +  R2  +  R3  (82) 
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o  -  {-{Yo*3(0)  +  «oK*4(0)> 

+  qi(1)  -  +  P12(1)*4(0)  +  P13(1)Y3(1)  +  P14(1)Y4<1))  }6Y3<°) 

+  (q2^^  “  [p21^^3^^  +  P22^^4^^  +  P23^3^^  +  P24^^4^^1  JdY^®^ 

K-l 

+  I  Uqi(mH)  +  q3Cm>  ] 

m»l 

-  [pu(nH-l)Y3(m)  +  p12(m+l)Y4(m)  +  pu(n»H)Y3(®H)  +  PiA(e*1) 

-  lP3l<m>Y3<m-1>  +  p32(m)Y4(m_1>  +  P33(m)Y3(m>  +  P34(m)Y4(m>  1  }<5Y3<ni> 

K-l 

+  I  {q2(nrfl)  +  q4<m>  1 

m-1 

-  [p2i^Bri‘1H3^m)  +  P22^®4‘1^4^m)  +  p23(nrt-l)Y3<nrt-l)  +  p24^ari‘1^Y4^ni+1^  I 

-  lP41<m)Y3<®-l>  +  P42(m)Y4(m_1>  +  P43(n,)*3(m)  +  P44<m>Y4<®>)  }6Y4<m>  (83) 


RECURSIVE  SOLUTIONS 

Since  the  variations  6Y3(0) ,  6Y4^)  ,  5Y3^m^  »  and  fiY^®)  in  Eq.  (83)  are 
all  arbitrary,  the  coefficients  of  all  these  variations  must  vanish.  Me  first 
take  the  coefficients  of  the  variations  6Y3^)  and  6Y4'®). 


P13^  Pl4^ 

y3U) 

P23^L)  P24^ 

Y4^> 

_  _ 

(_p11(l)_Yo)(-pi2(1)-aoK) 

y3(°) 

+ 

Vl) 

(_P2i(l)  (“P22^ 

Y4<°> 

q2(U 

_ 

—  _ 

(84) 

It  is  noted  that  Y3^  and  Y4^)  are  the  initial  conditions  of  the  problem 
that  is  from  Eq.  (67)  and  (46). 

Y3<°>  -  y0  (85) 


Y4<°> 


Yo 


)  1  dY 

-  K~lyQf  -  R'1  — 


dt 


(86) 
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We  can  solve  for  Y3O)  and  Y4(l)  in  terms  of  these  initial  conditions  by 


inverting  the  two  by  two  matrix  in  Eq.  (84). 

_  _  _  _  /"_  _  _  _  _  _ 

Y3(D  Pi3(1)  Pi4(1)  J  (-Pn(1)-Y0)  (-P12“«oK)  y0  qi(1) 

■  <  + 

Y/. ( 1)  p23^  J  P24^ ^  I  (  — o>i‘>^^)  K~^v'^  O')^  ^ 


For  a  general  case  where  m  =  1,2,...K-1,  we  have  by  setting  the 
coefficients  of  <5Y3^m)  and  6Y4^m^  in  Eq.  (83)  to  zero. 


Y3(ari-1) 

Pl3(mfl)  Pl4(nri_1) 

Y4(nri-D 

P23(mfl)  P24(mH) 

<Pll(n,+1)  +  P33(m))  (Pl2(lttfl)  +  P34(m>)  Y3(») 

(P21(n,+1)  +  P43U))  (P22(“ri'1)  +  P24U)>  *4(n) 


P3l(m)  P32^m^  Y3^n~1)  q^(nrt-l) 

+ 

P41(al)  P42U)  Y4<*“1)  q2^nH'l> 


We  solve  the  above  equation  recursively  for  Y3^nr*'^  and  Y4^m+'^.  Starting 
with  s  *  1  we  have  the  initial  conditions  Y3W)  and  Y^O)  and  the  solutions 
from  Eq.  (87)  for  Y3^  and  Y4^) .  These  values  are  substituted  into  Eq.  (88) 
to  obtain  ?3^)  and  Y4^) .  pr0m  the  values  of  Y3^,  Y4^\  Y3^)  f  and  Y4^) 
one  can  determine  Y3^)  and  Y4^).  This  procedure  continues  until  we  obtain 
Y3^K)  and  Y4^)  which  are  the  final  values  of  the  problem. 
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NUMERICAL  RESULTS  AND  DISCUSSION 

The  analysis  presented  In  previous  sections  will  now  be  tested  by  way  of 
some  numerical  examples.  Let  us  consider  a  simple  oscillator  subjected  to  a 
harmonic  force.  The  differential  equation  can  be  written  as 

ray  +  ky  »  f0  cos  u>ft  0  <  t  <  T  (89a) 

where  T  is  some  finite  time  of  interest  and  a  dot  ( *)  denotes  differentiation 
with  respect  to  time.  The  initial  conditions  are 

y(0)  -  y0  and  y(0)  -  y0  (89b) 

The  system  of  Eqs.  (89a)  and  (89b)  is  normalized  with  respect  to  T  and  it 
becomes 

y*"  +  k*y*  «  f*  cos  uf*t*  0  <  t*  <  1  (90a) 

and 

y*(0)  *  yo*  and  y*(0)  »•  yo*  (90b) 

Through  the  following  change  of  parameters 

t  dt 

t*  -  -  ,  dt*  -  — 

T  T 

y* ( t* )  *  y(t)  ,  y*'(t*)  »  T  —  (91) 

dt 

k*  -  kTz/m  ,  f*  -  f0T2/m  ,  uif*  **  ufT 

yo*  •  yo  »  yi*  *  Tyi 


Comparing  Eq.  (90a)  with  Eqs.  (24)  and  (1),  one  has 

a  “  constant  -  1  ,  e  *  -1 

8*0  ,  y  ■  0  and  Q  »  -f*  cos  a>f*t* 


(92) 
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From  the  data  presented  here,  we  further  set 


m  -  i  .0  ,  k  *»  1 .0  ,  t0  »  1.0  ,  uif  -  0.5 

The  parameter  T  is  given  for  each  set  of  sample  calculations. 

First,  Eq.  (84)  can  be  used  exclusively  to  obtain  all  the  solutions. 

This  is  demonstrated  in  Tables  I  through  III.  In  these  tables  T  has  taken  to 
be  ten,  five,  and  two,  respectively  and  the  number  of  steps  for  all  cases  is 
taken  to  be  ten.  Both  y(t)  and  y(t)  are  shown  and  the  exact  solutions  are 

given  in  parentheses  for  comparison.  It  is  clear  that  the  results  are 

convergent,  i.e.,  they  are  improved  as  the  interval  of  time  is  decreased. 

TABLE  I.  SOLUTION  TO  A  FORCED  VIBRATION  PROBLEM  OF  A  SIMPLE  OSCILLATOR 

(0  <  t  <  10,  Ten  Steps.  Exact  Solution  Shown  in  Parenthesis) 


rt 

y(t) 

y(t) 

0 

1.0000 

(Given) 

1.000 

(Given) 

2.0 

1.7590 

(  1.7684) 

-0.711 

(-0.674) 

4.0 

-1.1495 

(-1.0938) 

-1.450 

(-1.512) 

6.0 

-1.8534 

(-1.9195) 

0.867 

(  0.773) 

8.0 

0.2261 

(  0.1663) 

0.564 

(  0.689) 

10.0 

-0.0531 

(  0.1139) 

-0.404 

(-0.381) 
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TABLE  II.  SOLUTIONS  TO  A  FORCED  VIBRATION  PROBLEM  OF  A  SIMPLE  OSCILLATOR 


(0  <  t  <  3,  Ten  Steps.  Exact  Solutions  Shovm  in  Parentnesis) 


t 

y(t) 

y(t) 

0 

1 .0000 

(Given) 

1 .0000 

(Given) 

1.0 

1.8314 

(  1.8315) 

0.4991 

(  .5012) 

2.0 

1.7646 

(  1.7684) 

-0.6828 

(-0.6740) 

3.0 

0.5536 

(  0.5654) 

-1.6161 

(-1.6079) 

4.0 

-1.1074 

(-1.09^8) 

-1.5060 

(-1.5121) 

EH 

-2.1221 

(-2.1217) 

-0.4129 

(-0.4350) 

TABLE  III.  SOLUTIONS  TO  A  FORCED  VIBRATION  PROBLEM  OF  A  SIMPLE  OSCILLATOR 
(0  <  t  <  2.0,  Ten  Steps.  Exact  Solutions  Shown  in  Parenthesis) 


t 

y(t) 

y(t) 

g 

1 .0000 

(Given) 

1 .0000 

(Given) 

1.3892 

(1.3892) 

0.9184 

(  .9184) 

0.8 

1.7132 

(1.7132) 

0.6760 

(-0.6760) 

1.2 

1.9116 

(1.9117) 

0.2961 

(  0.2966) 

1.6 

1.9379 

(1.9382) 

-0.1752 

(-0.1742) 

2.0 

1.7676 

(1.7684) 

-0.6754 

(-0.6740) 

21 


Some  discussion  on  che  present  formulation  compared  with  previous  work3,4 


is  in  order  here.  In  previous  work  on  unconstrained,  adjoint  variational 
formulation,  the  point  of  emphasis  was  to  free  the  requirements  of  satisfying 
any  of  the  initial  conditions  and  to  let  the  approximate  solution  converge  to 
them.  In  the  present  analysis  it  is  shown  that  the  far  end  conditions  need 
not  be  considered  in  a  variational  formulation  of  approximate  solutions.  A 
more  detailed  comparison  in  terms  of  numerical  convergence,  competency, 
efficiency,  etc.  is  planned. 


3j.  J.  Wu,  ’’Solutions  to  Initial  Value  Problems  By  Use  of  Finite-Element- 
Unconstrained  Variational  Formulations,”  Journal  of  Sound  and  Vibration,  53 , 
1977,  pp.  344-356. 

^J.  J.  Wu  and  T.  E.  Slmklns,  "A  Numerical  Comparison  Between  Two  Unconstrained 
Variational  Formulations,”  Journal  of  Sound  and  Vibration,  72^,  1980,  pp. 
491-506. 
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